# Paper: Is it still the economy? Economic voting in polarized politics
# County-level simulation
# Replication: Figure 5
remove(list = ls())

library(conflicted)
library(tidyverse)
conflict_prefer("filter","dplyr")
library(ggplot2)
library(dplyr)
library(here)
conflict_prefer("here", "here")
library(systemfit)

load(here("Final_Paper", "R&R", "CountyLevel.RData"))

lia1 <- inc_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) + 
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff +
  inc_abs_lag + state + rep_inc + incumbent_cand
lia2 <- opp_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) +
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff +
  opp_abs_lag + state + rep_inc + incumbent_cand
Systemia <- list(lia1, lia2)
modelia <- systemfit(Systemia, method = "SUR", data = model_data)
summary(modelia)

lio1 <- inc_opp ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) + 
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff + 
  inc_opp_lag + state + rep_inc + incumbent_cand
lio2 <- abs_opp ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) + 
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff + 
  abs_opp_lag + state + rep_inc + incumbent_cand
Systemio <- list(lio1, lio2)
modelio <- systemfit(Systemio, method = "SUR", data = model_data)
summary(modelio)



################################################################################
################################################################################
################################################################################
###################################Bootstrap#################################### 
################################################################################
################################################################################
################################################################################

glimpse(model_data)
summary(model_data$unemp_rate)
sd(model_data$unemp_rate) #2.679395
sd(model_data$change_unemprate) #0.3491287
table(model_data$state)
sample(model_data$state, 1)
table(model_data$elitepolar_house)
summary(model_data$elitepolar_house)

# within state variation in unemployment rate
wstate_unemp <- lm(unemp_rate ~ state, model_data)
summary(wstate_unemp)
model_data <- model_data %>% mutate(within_unemp = resid(wstate_unemp))
sd(model_data$within_unemp)#2.357402
sd(model_data$unemp_rate)  #2.679395
mean(model_data$unemp_rate)#6.121542

model_data2 <- model_data %>% 
  select(year, state, unemp_rate, elitepolar_house, defl_pcincome, 
         share_black, share_hisp, state_partiesdiff,
         share_asian, share_young, share_elder, 
         tot_pop, share_college, share_manuf,
         rural_urban, inc_abs, inc_abs_lag, opp_abs, 
         opp_abs_lag, rep_inc, inc_opp, incumbent_cand,
         inc_opp_lag, abs_opp, abs_opp_lag, within_unemp)
#model_data1 <- model_data %>% filter(!(year == 1992))
glimpse(model_data2)



rep <- length(model_data2$inc_abs)
indices_rep <- 1:rep
#results_rep <- as.data.frame(matrix(data = NA, nrow = 134, ncol = 100))
results_absr <- as.data.frame(matrix(data = NA, nrow = 146, ncol = 1000))
results_transf_unemp <- as.data.frame(matrix(data = NA, nrow = 1000, ncol = 6))
#results_transf_pol <- as.data.frame(matrix(data = NA, nrow = 100, ncol = 3))
#results_transf_int <- as.data.frame(matrix(data = NA, nrow = 100, ncol = 3))
#sample_indices <- sample(indices_rep, rep, replace = T)
#bootstrap_data <- rep_data[sample_indices, ]

#abstention as the baseline
summary(modelia)
summary(modelia)$coefficients[2,  ] 
#unemp_rate                  -0.094507080  0.006168638 -15.32057 < 2.22e-16 ***
summary(modelia)$coefficients[73,  ] 
#unemp_rate:elitepolar_house  0.087366390  0.007895326  11.06558 < 2.22e-16 ***
table(model_data2$rural_urban) #1st eq.: rural_urban == 7
summary(modelia)$coefficients[18,  ] 
#as.factor(rural_urban)7      0.071772099  0.015980645   4.49119 7.1147e-06 ***
summary(modelia)$coefficients[58,  ] #1st eq.: PA state
#statePA                     -0.015799123  0.016446643  -0.96063  0.3367481 
summary(modelia)$coefficients[91,  ]  #2nd eq: rural_urban == 7
#as.factor(rural_urban)7      0.01440144  0.01669160   0.86280 0.38825817
summary(modelia)$coefficients[131,  ] #2nd eq: PA state
#statePA                      0.03946169  0.01717672   2.29739 0.02160467 *


#randomizer <- c(1:1000000)
#sample(randomizer, 1)
#remove(bootstrap_data, sample_indices)
set.seed(251827)
for(i in 1:1000) {
  
  # Resampling w/ replacement
  tryCatch({sample_indices <- sample(indices_rep, rep, replace = T)
  bootstrap_data <- model_data2[sample_indices, ]
  
  #SUR Model 1: Abstention as the baseline category
  modelia <- systemfit(list(inc_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
                              log(defl_pcincome) +
                              share_black + share_hisp + share_asian + share_young + 
                              share_elder + log(tot_pop) + 
                              log1p(share_college) + log1p(share_manuf) +
                              as.factor(rural_urban) + state_partiesdiff +
                              inc_abs_lag + state + rep_inc + incumbent_cand,
                            opp_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
                              log(defl_pcincome) +
                              share_black + share_hisp + share_asian + share_young + 
                              share_elder + log(tot_pop) +
                              log1p(share_college) + log1p(share_manuf) +
                              as.factor(rural_urban) + state_partiesdiff +
                              opp_abs_lag + state + rep_inc + incumbent_cand), 
                       method = "SUR", data = bootstrap_data)
  
  # storing coefficients
  results_absr[, i] <- as.data.frame(summary(modelia)$coefficients[,  1])
  names(results_absr)[i] <- paste("coeff", i, sep = "_")
  
  #results_rep[, i] <- as.data.frame(summary(modelr)$coefficients[,  1])
  #names(results_rep)[i] <- paste("coeff", i, sep = "_")
  
  ####### 1 SD shock in unemployment with low polarization
  denom1 <- 1 + exp(modelia$coefficients[1] + 
                      modelia$coefficients[2]*(mean(model_data2$unemp_rate) + 
                                                 sd(model_data2$within_unemp)) + 
                      modelia$coefficients[3]*0.60 + #polarization
                      modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
                      modelia$coefficients[5]*mean(model_data2$share_black) +
                      modelia$coefficients[6]*mean(model_data2$share_hisp) +
                      modelia$coefficients[7]*mean(model_data2$share_asian) +
                      modelia$coefficients[8]*mean(model_data2$share_young) +
                      modelia$coefficients[9]*mean(model_data2$share_elder) +
                      modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
                      modelia$coefficients[11]*log(mean(model_data2$share_college)) +
                      modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
                      modelia$coefficients[19]*1 + #rural_urban == 7
                      modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
                      modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
                      modelia$coefficients[59]*1 + #Pennsylvania
                      modelia$coefficients[72]*1 +
                      modelia$coefficients[73]*1 +
                      modelia$coefficients[74]*(mean(model_data2$unemp_rate) + 
                                                  sd(model_data2$within_unemp))*0.60) +
    exp(modelia$coefficients[75] + 
          modelia$coefficients[76]*(mean(model_data2$unemp_rate) + 
                                      sd(model_data2$within_unemp))  + 
          modelia$coefficients[77]*0.60 + 
          modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
          modelia$coefficients[79]*mean(model_data2$share_black) +
          modelia$coefficients[80]*mean(model_data2$share_hisp) +
          modelia$coefficients[81]*mean(model_data2$share_asian) +
          modelia$coefficients[82]*mean(model_data2$share_young) +
          modelia$coefficients[83]*mean(model_data2$share_elder) +
          modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
          modelia$coefficients[85]*log(mean(model_data2$share_college)) +
          modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
          modelia$coefficients[93]*1 + #rural_urban == 7
          modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
          modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
          modelia$coefficients[124]*1 +#Pennsylvania
          modelia$coefficients[146]*1 +
          modelia$coefficients[147]*1 +
          modelia$coefficients[148]*(mean(model_data2$unemp_rate)+ 
                                       sd(model_data2$within_unemp))*0.60)
  
  denom2 <- 1 + exp(modelia$coefficients[1] + 
                      modelia$coefficients[2]*mean(model_data2$unemp_rate) + 
                      modelia$coefficients[3]*0.60 + #polarization
                      modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
                      modelia$coefficients[5]*mean(model_data2$share_black) +
                      modelia$coefficients[6]*mean(model_data2$share_hisp) +
                      modelia$coefficients[7]*mean(model_data2$share_asian) +
                      modelia$coefficients[8]*mean(model_data2$share_young) +
                      modelia$coefficients[9]*mean(model_data2$share_elder) +
                      modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
                      modelia$coefficients[11]*log(mean(model_data2$share_college)) +
                      modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
                      modelia$coefficients[19]*1 + #rural_urban == 7
                      modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
                      modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
                      modelia$coefficients[59]*1 + #Pennsylvania
                      modelia$coefficients[72]*1 +
                      modelia$coefficients[73]*1 +
                      modelia$coefficients[74]*mean(model_data2$unemp_rate)*0.60) +
    exp(modelia$coefficients[75] + 
          modelia$coefficients[76]*mean(model_data2$unemp_rate)  + 
          modelia$coefficients[77]*0.60 + 
          modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
          modelia$coefficients[79]*mean(model_data2$share_black) +
          modelia$coefficients[80]*mean(model_data2$share_hisp) +
          modelia$coefficients[81]*mean(model_data2$share_asian) +
          modelia$coefficients[82]*mean(model_data2$share_young) +
          modelia$coefficients[83]*mean(model_data2$share_elder) +
          modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
          modelia$coefficients[85]*log(mean(model_data2$share_college)) +
          modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
          modelia$coefficients[93]*1 + #rural_urban == 7
          modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
          modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
          modelia$coefficients[124]*1 +#Pennsylvania
          modelia$coefficients[146]*1 +
          modelia$coefficients[147]*1 +
          modelia$coefficients[148]*mean(model_data2$unemp_rate)*0.60)
  
  
  #Change in Incumbent vote share
  results_transf_unemp[i, 1] <- (exp(modelia$coefficients[1] + 
                                       modelia$coefficients[2]*(mean(model_data2$unemp_rate) + 
                                                                  sd(model_data2$within_unemp)) + 
                                       modelia$coefficients[3]*0.60 + #polarization
                                       modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
                                       modelia$coefficients[5]*mean(model_data2$share_black) +
                                       modelia$coefficients[6]*mean(model_data2$share_hisp) +
                                       modelia$coefficients[7]*mean(model_data2$share_asian) +
                                       modelia$coefficients[8]*mean(model_data2$share_young) +
                                       modelia$coefficients[9]*mean(model_data2$share_elder) +
                                       modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
                                       modelia$coefficients[11]*log(mean(model_data2$share_college)) +
                                       modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
                                       modelia$coefficients[19]*1 + #rural_urban == 7
                                       modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
                                       modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
                                       modelia$coefficients[59]*1 + #Pennsylvania
                                       modelia$coefficients[72]*1 +
                                       modelia$coefficients[73]*1 +
                                       modelia$coefficients[74]*(mean(model_data2$unemp_rate) + 
                                                                   sd(model_data2$within_unemp))*0.60) / denom1) -
    (exp(modelia$coefficients[1] + 
           modelia$coefficients[2]*mean(model_data2$unemp_rate) + 
           modelia$coefficients[3]*0.60 + #polarization
           modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
           modelia$coefficients[5]*mean(model_data2$share_black) +
           modelia$coefficients[6]*mean(model_data2$share_hisp) +
           modelia$coefficients[7]*mean(model_data2$share_asian) +
           modelia$coefficients[8]*mean(model_data2$share_young) +
           modelia$coefficients[9]*mean(model_data2$share_elder) +
           modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
           modelia$coefficients[11]*log(mean(model_data2$share_college)) +
           modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
           modelia$coefficients[19]*1 + #rural_urban == 7
           modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
           modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
           modelia$coefficients[59]*1 + #Pennsylvania
           modelia$coefficients[72]*1 +
           modelia$coefficients[73]*1 +
           modelia$coefficients[74]*mean(model_data2$unemp_rate)*0.60) / denom2) 
  
  #Change in Opposition vote share
  results_transf_unemp[i, 2]  <- (exp(modelia$coefficients[75] + 
                                        modelia$coefficients[76]*(mean(model_data2$unemp_rate) + 
                                                                    sd(model_data2$within_unemp))  + 
                                        modelia$coefficients[77]*0.60 + 
                                        modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
                                        modelia$coefficients[79]*mean(model_data2$share_black) +
                                        modelia$coefficients[80]*mean(model_data2$share_hisp) +
                                        modelia$coefficients[81]*mean(model_data2$share_asian) +
                                        modelia$coefficients[82]*mean(model_data2$share_young) +
                                        modelia$coefficients[83]*mean(model_data2$share_elder) +
                                        modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
                                        modelia$coefficients[85]*log(mean(model_data2$share_college)) +
                                        modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
                                        modelia$coefficients[93]*1 + #rural_urban == 7
                                        modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
                                        modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
                                        modelia$coefficients[124]*1 +#Pennsylvania
                                        modelia$coefficients[146]*1 +
                                        modelia$coefficients[147]*1 +
                                        modelia$coefficients[148]*(mean(model_data2$unemp_rate)+ 
                                                                     sd(model_data2$within_unemp))*0.60) / denom1) - 
    (exp(modelia$coefficients[75] + 
           modelia$coefficients[76]*mean(model_data2$unemp_rate)  + 
           modelia$coefficients[77]*0.60 + 
           modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
           modelia$coefficients[79]*mean(model_data2$share_black) +
           modelia$coefficients[80]*mean(model_data2$share_hisp) +
           modelia$coefficients[81]*mean(model_data2$share_asian) +
           modelia$coefficients[82]*mean(model_data2$share_young) +
           modelia$coefficients[83]*mean(model_data2$share_elder) +
           modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
           modelia$coefficients[85]*log(mean(model_data2$share_college)) +
           modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
           modelia$coefficients[93]*1 + #rural_urban == 7
           modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
           modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
           modelia$coefficients[124]*1 +#Pennsylvania
           modelia$coefficients[146]*1 +
           modelia$coefficients[147]*1 +
           modelia$coefficients[148]*mean(model_data2$unemp_rate)*0.60) / denom2)
  
  #Change in proportion of abstention
  results_transf_unemp[i, 3]  <- (1/denom1) -
    (1/denom2)
  
  ####### 1 SD shock in unemployment with high polarization
  denom3 <- 1 + exp(modelia$coefficients[1] + 
                      modelia$coefficients[2]*(mean(model_data2$unemp_rate) + 
                                                 sd(model_data2$within_unemp)) + 
                      modelia$coefficients[3]*0.85 + #polarization
                      modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
                      modelia$coefficients[5]*mean(model_data2$share_black) +
                      modelia$coefficients[6]*mean(model_data2$share_hisp) +
                      modelia$coefficients[7]*mean(model_data2$share_asian) +
                      modelia$coefficients[8]*mean(model_data2$share_young) +
                      modelia$coefficients[9]*mean(model_data2$share_elder) +
                      modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
                      modelia$coefficients[11]*log(mean(model_data2$share_college)) +
                      modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
                      modelia$coefficients[19]*1 + #rural_urban == 7
                      modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
                      modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
                      modelia$coefficients[59]*1 + #Pennsylvania
                      modelia$coefficients[72]*1 +
                      modelia$coefficients[73]*1 +
                      modelia$coefficients[74]*(mean(model_data2$unemp_rate) + 
                                                  sd(model_data2$within_unemp))*0.85) +
    exp(modelia$coefficients[75] + 
          modelia$coefficients[76]*(mean(model_data2$unemp_rate) + 
                                      sd(model_data2$within_unemp))  + 
          modelia$coefficients[77]*0.85 + 
          modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
          modelia$coefficients[79]*mean(model_data2$share_black) +
          modelia$coefficients[80]*mean(model_data2$share_hisp) +
          modelia$coefficients[81]*mean(model_data2$share_asian) +
          modelia$coefficients[82]*mean(model_data2$share_young) +
          modelia$coefficients[83]*mean(model_data2$share_elder) +
          modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
          modelia$coefficients[85]*log(mean(model_data2$share_college)) +
          modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
          modelia$coefficients[93]*1 + #rural_urban == 7
          modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
          modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
          modelia$coefficients[124]*1 +#Pennsylvania
          modelia$coefficients[146]*1 +
          modelia$coefficients[147]*1 +
          modelia$coefficients[148]*(mean(model_data2$unemp_rate)+ 
                                       sd(model_data2$within_unemp))*0.85)
  
  denom4 <- 1 + exp(modelia$coefficients[1] + 
                      modelia$coefficients[2]*mean(model_data2$unemp_rate) + 
                      modelia$coefficients[3]*0.85 + #polarization
                      modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
                      modelia$coefficients[5]*mean(model_data2$share_black) +
                      modelia$coefficients[6]*mean(model_data2$share_hisp) +
                      modelia$coefficients[7]*mean(model_data2$share_asian) +
                      modelia$coefficients[8]*mean(model_data2$share_young) +
                      modelia$coefficients[9]*mean(model_data2$share_elder) +
                      modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
                      modelia$coefficients[11]*log(mean(model_data2$share_college)) +
                      modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
                      modelia$coefficients[19]*1 + #rural_urban == 7
                      modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
                      modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
                      modelia$coefficients[59]*1 + #Pennsylvania
                      modelia$coefficients[72]*1 +
                      modelia$coefficients[73]*1 +
                      modelia$coefficients[74]*mean(model_data2$unemp_rate)*0.85) +
    exp(modelia$coefficients[75] + 
          modelia$coefficients[76]*mean(model_data2$unemp_rate)  + 
          modelia$coefficients[77]*0.85 + 
          modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
          modelia$coefficients[79]*mean(model_data2$share_black) +
          modelia$coefficients[80]*mean(model_data2$share_hisp) +
          modelia$coefficients[81]*mean(model_data2$share_asian) +
          modelia$coefficients[82]*mean(model_data2$share_young) +
          modelia$coefficients[83]*mean(model_data2$share_elder) +
          modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
          modelia$coefficients[85]*log(mean(model_data2$share_college)) +
          modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
          modelia$coefficients[93]*1 + #rural_urban == 7
          modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
          modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
          modelia$coefficients[124]*1 +#Pennsylvania
          modelia$coefficients[146]*1 +
          modelia$coefficients[147]*1 +
          modelia$coefficients[148]*mean(model_data2$unemp_rate)*0.85)
  
  
  #Change in Incumbent vote share
  results_transf_unemp[i, 4] <- (exp(modelia$coefficients[1] + 
                                       modelia$coefficients[2]*(mean(model_data2$unemp_rate) + 
                                                                  sd(model_data2$within_unemp)) + 
                                       modelia$coefficients[3]*0.85 + #polarization
                                       modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
                                       modelia$coefficients[5]*mean(model_data2$share_black) +
                                       modelia$coefficients[6]*mean(model_data2$share_hisp) +
                                       modelia$coefficients[7]*mean(model_data2$share_asian) +
                                       modelia$coefficients[8]*mean(model_data2$share_young) +
                                       modelia$coefficients[9]*mean(model_data2$share_elder) +
                                       modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
                                       modelia$coefficients[11]*log(mean(model_data2$share_college)) +
                                       modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
                                       modelia$coefficients[19]*1 + #rural_urban == 7
                                       modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
                                       modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
                                       modelia$coefficients[59]*1 + #Pennsylvania
                                       modelia$coefficients[72]*1 +
                                       modelia$coefficients[73]*1 +
                                       modelia$coefficients[74]*(mean(model_data2$unemp_rate) + 
                                                                   sd(model_data2$within_unemp))*0.85) / denom3) -
    (exp(modelia$coefficients[1] + 
           modelia$coefficients[2]*mean(model_data2$unemp_rate) + 
           modelia$coefficients[3]*0.85 + #polarization
           modelia$coefficients[4]*log(mean(model_data2$defl_pcincome)) + 
           modelia$coefficients[5]*mean(model_data2$share_black) +
           modelia$coefficients[6]*mean(model_data2$share_hisp) +
           modelia$coefficients[7]*mean(model_data2$share_asian) +
           modelia$coefficients[8]*mean(model_data2$share_young) +
           modelia$coefficients[9]*mean(model_data2$share_elder) +
           modelia$coefficients[10]*log(mean(model_data2$tot_pop)) +
           modelia$coefficients[11]*log(mean(model_data2$share_college)) +
           modelia$coefficients[12]*log1p(mean(model_data2$share_manuf)) +
           modelia$coefficients[19]*1 + #rural_urban == 7
           modelia$coefficients[22]*mean(model_data2$state_partiesdiff, na.rm = T) +
           modelia$coefficients[23]*mean(model_data2$inc_abs_lag, na.rm = T) +
           modelia$coefficients[59]*1 + #Pennsylvania
           modelia$coefficients[72]*1 +
           modelia$coefficients[73]*1 +
           modelia$coefficients[74]*mean(model_data2$unemp_rate)*0.85) / denom4) 
  
  #Change in Opposition vote share
  results_transf_unemp[i, 5]  <- (exp(modelia$coefficients[75] + 
                                        modelia$coefficients[76]*(mean(model_data2$unemp_rate) + 
                                                                    sd(model_data2$within_unemp))  + 
                                        modelia$coefficients[77]*0.85 + 
                                        modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
                                        modelia$coefficients[79]*mean(model_data2$share_black) +
                                        modelia$coefficients[80]*mean(model_data2$share_hisp) +
                                        modelia$coefficients[81]*mean(model_data2$share_asian) +
                                        modelia$coefficients[82]*mean(model_data2$share_young) +
                                        modelia$coefficients[83]*mean(model_data2$share_elder) +
                                        modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
                                        modelia$coefficients[85]*log(mean(model_data2$share_college)) +
                                        modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
                                        modelia$coefficients[93]*1 + #rural_urban == 7
                                        modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
                                        modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
                                        modelia$coefficients[124]*1 +#Pennsylvania
                                        modelia$coefficients[146]*1 +
                                        modelia$coefficients[147]*1 +
                                        modelia$coefficients[148]*(mean(model_data2$unemp_rate)+ 
                                                                     sd(model_data2$within_unemp))*0.85) / denom3) - 
    (exp(modelia$coefficients[75] + 
           modelia$coefficients[76]*mean(model_data2$unemp_rate)  + 
           modelia$coefficients[77]*0.85 + 
           modelia$coefficients[78]*log(mean(model_data2$defl_pcincome)) + 
           modelia$coefficients[79]*mean(model_data2$share_black) +
           modelia$coefficients[80]*mean(model_data2$share_hisp) +
           modelia$coefficients[81]*mean(model_data2$share_asian) +
           modelia$coefficients[82]*mean(model_data2$share_young) +
           modelia$coefficients[83]*mean(model_data2$share_elder) +
           modelia$coefficients[84]*log(mean(model_data2$tot_pop)) +
           modelia$coefficients[85]*log(mean(model_data2$share_college)) +
           modelia$coefficients[86]*log1p(mean(model_data2$share_manuf)) +
           modelia$coefficients[93]*1 + #rural_urban == 7
           modelia$coefficients[96]*mean(model_data2$state_partiesdiff, na.rm = T) +
           modelia$coefficients[97]*mean(model_data2$opp_abs_lag, na.rm = T) +
           modelia$coefficients[124]*1 +#Pennsylvania
           modelia$coefficients[146]*1 +
           modelia$coefficients[147]*1 +
           modelia$coefficients[148]*mean(model_data2$unemp_rate)*0.85) / denom4)
  
  #Change in proportion of abstention
  results_transf_unemp[i, 6] <- (1/denom3) -
    (1/denom4)
  
  print(i)}, error=function(e){cat("Error in LU.dgC(a) : cs_lu(A) failed: near-singular A (or out of memory)",
                                   conditionMessage(e), "\n")}) 
}
#2 error

glimpse(results_transf_unemp)

############### Saving Coefficients
lia1 <- inc_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) + 
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff +
  inc_abs_lag + state + rep_inc + incumbent_cand
lia2 <- opp_abs ~ unemp_rate*elitepolar_house + #per_chg_wage*pol_median_dist1 + 
  log(defl_pcincome) +
  share_black + share_hisp + share_asian + share_young + share_elder + log(tot_pop) +
  log1p(share_college) + log1p(share_manuf) +
  as.factor(rural_urban) + state_partiesdiff +
  opp_abs_lag + state + rep_inc + incumbent_cand
Systemra <- list(lia1, lia2)
modelia <- systemfit(Systemra, method = "SUR", data = model_data)
summary(modelia)
coeffdata <- as.data.frame(cbind(names(modelia$coefficients), modelia$coefficients))
names(coeffdata) <- c("varnames", "mid")
coeffdata$mid <- as.numeric(coeffdata$mid)

glimpse(results_absr)
class(results_absr)
results_absr$varnames <- names(modelia$coefficients)
results_absr <- results_absr[ , colSums(is.na(results_absr)) == 0]
results_absr1 <- gather(results_absr, key = "coeff_", value = "test", coeff_1:coeff_1000)
glimpse(results_absr1)

results_absr1 <- results_absr1 %>%
  group_by(varnames) %>%
  summarise(sd = sd(test)) %>%
  ungroup()
glimpse(results_absr1)
glimpse(coeffdata)
results_absr1 <- left_join(results_absr1, coeffdata, by = "varnames")

results_absr1 <- results_absr1 %>% 
  mutate(equation = substr(varnames, 1, 3)) %>% 
  mutate(varnames = substr(varnames, 5, nchar(varnames))) %>% 
  mutate(lb = mid - 1.96*sd) %>% 
  mutate(ub = mid + 1.96*sd)

glimpse(results_absr1)
#write.csv(results_absr1, here("Final_Paper", "coeff_sim_unemp_pa.csv"))

############### Simulated vote share
glimpse(results_transf_unemp)

results_transf_unemp <- results_transf_unemp %>% 
  rename(inc_low = V1, opp_low = V2, abs_low = V3,
         inc_high = V4, opp_high = V5, abs_high = V6)

glimpse(results_transf_unemp)

results_transf_unemp <- results_transf_unemp[rowSums(is.na(results_transf_unemp)) == 0, ]
glimpse(results_transf_unemp)
#write.csv(results_transf_unemp, here("Final_Paper", "VoteShare_Sim_Unemp_pa.csv"))
#results_transf_unemp <- read.csv(here("Final_Paper", "VoteShare_Sim_Unemp.csv"))

party <- c("inc", "opp", "abs", "inc", "opp", "abs")
polar <- c(0.7, 0.7, 0.7, 0.8, 0.8, 0.8)
simul_data <- as.data.frame(matrix(NA, nrow = 6, ncol = 4))
simul_data <- simul_data %>% 
  rename(party = V1, polar = V2, mean = V3, sd = V4)
glimpse(simul_data)
simul_data[, 1] <- party
simul_data[, 2] <- polar
simul_data[1, 3] <- mean(results_transf_unemp$inc_low)
simul_data[2, 3] <- mean(results_transf_unemp$opp_low)
simul_data[3, 3] <- mean(results_transf_unemp$abs_low)
simul_data[4, 3] <- mean(results_transf_unemp$inc_high)
simul_data[5, 3] <- mean(results_transf_unemp$opp_high)
simul_data[6, 3] <- mean(results_transf_unemp$abs_high)
simul_data[1, 4] <- sd(results_transf_unemp$inc_low)
simul_data[2, 4] <- sd(results_transf_unemp$opp_low)
simul_data[3, 4] <- sd(results_transf_unemp$abs_low)
simul_data[4, 4] <- sd(results_transf_unemp$inc_high)
simul_data[5, 4] <- sd(results_transf_unemp$opp_high)
simul_data[6, 4] <- sd(results_transf_unemp$abs_high)

simul_data <- simul_data %>%
  mutate(lb = mean - 1.96*sd) %>% 
  mutate(ub = mean + 1.96*sd)

simul_data <- simul_data %>% 
  mutate(polar = as.factor(polar))

simul_data <- simul_data %>%
  mutate(party = ifelse(party == "inc", "Incumbent", party),
         party = ifelse(party == "opp", "Opposition", party),
         party = ifelse(party == "abs", "Abstention", party))

glimpse(simul_data)
fig_incu <- simul_data %>% 
  mutate(party = factor(party, levels = c("Incumbent",
                                          "Opposition",
                                          "Abstention"))) %>% 
  ggplot(aes(shape = polar, color = polar, fill = polar)) +
  geom_hline(yintercept = 0, colour = "gray37", lty = 2, size = 1) +
  geom_pointrange(aes(x = party, y = mean, ymin = lb,
                      ymax = ub),
                  lwd = 0.8, position = position_dodge(width = .5)) + 
  scale_colour_manual(name = "Polarization:",
                      values = c("#009999", "#990000"), 
                      labels = c("0.60", "0.85")) +
  scale_shape_manual(name = "Polarization:",
                     labels = c("0.60", "0.85"),
                     values = c(24, 25)) +
  scale_fill_manual(name = "Polarization:",
                    labels = c("0.60", "0.85"),
                    values = c("#009999", "#990000")) +
  theme_bw() + ylim(-0.025, 0.02) +
  labs(y = "Difference in Predicted Share of Voters") +#, 
  #   caption = "Effect of +1 Standard Deviation Change in Unemployment Rate") +
  xlab(NULL) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 13),
        legend.position = "bottom",
        legend.direction = "horizontal",
        title = element_text(size = 14))
fig_incu
